%% report Field that must contain the following 
% 
% 3) *.stateSpectrum*  Spectrum 
% 
% 2) *.stateIRF* Impulse responses 
% 
% 1) *.stateMom* Will be usded for all moments (stds and cross-correlations)
% including variance decompositions, and also Historical Decompositions 
%  
% All of the above will be computed automatically for the observables 
% so no need to define them. 
% 
% =========================================================================
cucd=changeRoot(cucd,'exog'); 
outpath=changeRoot(outpath,'exog');
loadpath=changeRoot(loadpath,'exog');



%% 2. Trim unfinished iterations 
posTrim=isinf(abs(logpostd_mat(:,2)));
if ~isempty(posTrim);
    logpostd_mat=logpostd_mat(posTrim==0,:);
    loglikel_mat=loglikel_mat(posTrim==0,:);
    logprior_mat=logprior_mat(posTrim==0,:);
    parmodemat=parmodemat(:,posTrim==0);
    stdLoopMat=stdLoopMat(:,posTrim==0);
    optimat=optimat(:,posTrim==0);  
    hessianMat=hessianMat(:,:,posTrim==0); 
end
clear posTrim;
    
%% 3. *orderDescend* is the indicator of the ordered modes 
cd(cucd);
[~,orderDescend]=sort(logpostd_mat(:, 2),'descend');
parModeMat=parmodemat(:, orderDescend);
parMode=parModeMat(:,1);
logLikelMat=loglikel_mat(orderDescend,:);
logPostMat=logpostd_mat(orderDescend,:);
hessianMat=hessianMat(:,:,orderDescend);
if logPostMat(1,2) < max(logPostMat(:,2) )
    error('Wrong Sorting of Results')
end
logPriorMat=logprior_mat(orderDescend,:);
stdLoopMat =stdLoopMat(:,orderDescend);
optimMat=optimat(:,orderDescend);
%clear orderDescend; %% Neded for Tables below, delete later

%% 4. Store to be used with zb_historical 
save2Priorstrv(parModeMat,Y,sample,Ynames,parnames,funcmod,addsol,solveopt,outpath,sampleSS,flagSplit); 

%% 6. Compute Hessian 
if flagStru.hessianDone==0
    if isfield(flagStru,'noHessian')==true && flagStru.noHessian==1
        stdHFD=nan( length(parposest) , 1 );
    else        
        disp('Begin Finite Difference Hessian');
        estimStructure.parVec(parposest)=zeros(length(parposest),1);
        estimStructure.LPosterior=logPostMat(1,2);
        [HFD,stdHFD,flagNotPD]=feval(@estimaHessian,parMode(parposest),Y,funcmod,solveopt,addsol,prpar,prss,estimStructure);
        disp('End Finite Difference Hessian');
        crtcell([parnames(parposest) num2cprec([stdHFD(:) stdLoopMat(:,1)]) ]);
        flagStru.hessianDone=1;
        close all;
        cd(outpath);
        save workspace;
        cd(cucd);
        disp('Hessian Done');
    end 
end

% ========================================================================
%% 7. Create ALL OUTPUT CELLS 
if flagStru.tableDone==0;
    estimaOutTabs;
    flagStru.tableDone=1;
end

%% 5. If ADD2FUNC exists, change function handle to this model after
if exist('add2Func','var')==true
if isempty(add2Func)==false && flagStru.handleChanged==0; 
    funcmod=handleChange(add2Func,funcmod,parMode,solveopt,addsol); 
    flagStru.handleChanged=1; 
end 
end


% ========================================================================
%% 8. Ploting STD, Vardecom, IRF, AC, KF state, KF Error for the best mode
if flag_mixfreq == 0;
    [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec,~, ssNames,stateNames,shockNamesStru]...
        =feval(funcmod,parMode,solveopt,addsol);
else
    disp('Re-do Completely!!!!'); 
    keyboard; 
    [GG, RR, CONS, eu, SDX, ZZ_zero, out, ssvec,~,ssNames,stateNames,shockNamesStru]...
        =feval(funcmod,parMode,solveopt,addsol);
    ZZ  =[ZZ_zero;out.AA];
    % added by XI Luo on 03/29/2011
    if ~isfield(out, 'case_tag') % if case_tag does not exist;
        % create one ourselves
        case_tag = nan(out.nzLF, 1);
        if out.flag_tagsum == 0;
            case_tag = 2*ones(out.nzLF, 1); % case 2: average
        elseif out.flag_tagsum == 1;
            for ii = 1:out.nzLF
                if max(out.AA(ii,:)) == 1/out.NUMIT;
                    case_tag(ii,1) = 2;
                elseif max(out.AA(ii,:)) == 1;
                    case_tag(ii,1) = 3;
                end
            end
        end
        out.case_tag = [-1*ones(out.nzHF, 1); case_tag]; % no transformation on high frequency
    end
end


%% ========================================================================
% Assignment of three key structures 
% 1. *nameStru:* Structure with the names 
%
% 2. report:* Series for which to report decomposition, IRFS and moments 
%
% 3. *momentStru:* Series for which to report moments 

%% *nameStru*
shockNames=shockNamesStru.long; 
[Znames,ZposInStates]=extractObsNames(ZZ(:,:,1),stateNames);
nameStru.stateNames=stateNames;
nameStru.shockNames=shockNames;
nameStru.obsNames=Znames;
%nameStru.reportObs=Znames;
%nameStru.reportStates=statesIRF;
%nameStru.NRep=1000; 
disp('Delete me later'); 
shockNames=shockNamesStru.long; 

%% *report*
[report,tempflag]=reportCheck(report,nameStru,sampleSS);  
flagStru=structureMerge(tempflag,flagStru,1); 
clear tempflag; 

[KFStru,logL]=feval(addsol.funcKfilter,parMode,funcmod,Y,trainvec,solveopt,addsol); 


%% Save all information of the model in structures 
estimaOutputStructures; 
disp('Saved model structures'); 


%% priorStru 
priorStru.distribution=prpar_dist; 
priorStru.mean=prpar_lcone; 
priorStru.std=prpar_lctwo; 
priorStru.lowerBound =prpar_lb; 
priorStru.upperBound =prpar_ub; 

%% postStru 
postStru.param=parMode; 
postStru.std=zeros(length(parMode),1); 
postStru.std(parposest)=stdLoopMat(:,1); 
postStru.logLikel=logLikelMat(1,2); 
postStru.logPost=logPostMat(1,2); 


namesTab.outputFile=['Report ',strdate]; 
namesTab.model=root; 
namesTab.specification=spec; 
namesTab.subfolder=subf; 
namesTab.param=parnames; 

if exist('description','var')==false 
    description=[]; 
end
cd(cucd); 
fclose('all'); 
cd(outpath);
fid=createReport(outpath,namesTab,priorStru,postStru,description); 
cd(cucd); 

%% names. 
% .outputFile: name of the file 
% .model: first layer 
% .specification: second layer 
% .subfolder: third layer 
% .param

% ========================================================================
%% 8. Spectrum  
cd(cucd); 
if flagStru.spectrumDone==0;
    if flagStru.stateSpectrum==false
        report.stateSpectrum=[];
    end        
    [spectrumCell,tableSpectrum]=estimaSpectrum(funcmod,parMode,solveopt,addsol,report.stateSpectrum,sampleSS,1000,outpath);
    clear formatStru;
    if length(tableSpectrum.colLabels) > 8
        formatStru.landscape=1;
    else
        formatStru.landscape=0;
    end
    fclose('all');
    cd(outpath)
    fid=fopen([namesTab.outputFile,'.tex'],'a+');    
    fprintf(fid,'\\pagebreak \n');
    latexTableCell(fid,num2cprec( tableSpectrum.dataMat(:,:,end) ,3 ),...
        tableSpectrum.rowLabels,tableSpectrum.colLabels,formatStru,{'Spectral Decomposition'},...
        'Series','Spectrum');
    clear tableSpectrum;
%     cd(outpath);
%     %xlswrite(xlsFilename,spectrumCell,'SpectralDecomposition');
%     cd(cucd);
     flagStru.spectrumDone=1;
    cd(outpath);
    save workspace;
    cd(cucd);
    disp('Spectrum Done');
end
fclose(fid); 


%% 
% ========================================================================
%% 8. Theoretical Moments 
% disp('Delete me later...'); 
% statesIRF=state_sel; 
estimaStates; 


%% Write to report
cd(outpath);
save workspace;
cd(cucd);




%% *momenStru* REMOVE, DONE IN momentsDraws
momentStru.NVecSim=[100 size(Y,1)]; 
momentStru.NRep=1000; 
momentStru.NVarDec=20; 
momentStru.NCorrel=4;
momentStru.NIRF=16; 
momentStru=momentsAssignDefaults(momentStru,report,nameStru); 



cd(cucd); 



%% First Sample Only  REMOVE, DONE IN momentsDraws
disp('Compile this into a single function'); 
momentOut=momentsDraws(GG(:,:,end),CONS(:,end),RR(:,:,end),SDX(:,:,end),ZZ(:,:,end),momentStru);
% USE momentOut from modelAnalysis

%% 8.a Impulse Responses REMOVE, DONE IN momentsDraws
nameStru.rowNames=report.stateIRF; 
nameStru.colNames=nameStru.shockNames; 
handlesIRF=plotIRF(momentOut.IRFStates,[],nameStru); 
compileGraphsPDF(handlesIRF,'IRF',outpath); 

%% 8.b Asymptotic Variance Decompositions Observables
tableStru.dataMat=zeros([size(momentOut.varDecHorObs(:,:,1)) 4]); 
tableStru.dataMat(:,:,1)=momentOut.varDecHorObs(:,:,1); 
tableStru.dataMat(:,:,2)=momentOut.varDecHorObs(:,:,4); 
tableStru.dataMat(:,:,3)=momentOut.varDecHorObs(:,:,8); 
tableStru.dataMat(:,:,4)=momentOut.varDecObs; 


%%

tableStru.rowLabels=Znames; 
tableStru.colLabels=shockNamesStru.short; 
tableStru.corner={'Observables'};
tableStru.tableHeader={'Variance Decomposition 1 Period Horizon',...
    'Variance Decomposition 2 Period Horizon',...
    'Variance Decomposition 4 Period Horizon',...
    'Variance Decomposition Asymptotic'};



if length(tableStru.colLabels) > 8
    formatStru.landscape=1;
end
fclose('all');
cd(outpath)
fid=fopen([namesTab.outputFile,'.tex'],'a+');
fprintf(fid,'\\pagebreak \n');
for counter=1:4
    latexTableCell(fid,num2cprec( tableStru.dataMat(:,:,counter) ,3 ),...
        tableStru.rowLabels, tableStru.colLabels,formatStru,tableStru.tableHeader(counter),...
        'Series',['Decomposition F Horizon=',num2str(counter)]);
    if counter==2 || counter ==4
        fprintf(fid,'\\pagebreak \n');
        
    end
end
fclose(fid); 

latexStru.printPDF=true;
latexStru.fileName=['Variance Decomposition Observables',strdate];
latexStru.outPath=outpath;
latexStru.fileExists=false;
latexStru.maxCols=8;
latexStru.maxRows=20;
latexStru.orientation='landscape';
multiPageTable(tableStru,latexStru);

%% 8.c Asymptotic Variance Decomposition States
tableStru.dataMat=zeros([size(momentOut.varDecHorStates(:,:,1)) 4]); 
tableStru.dataMat(:,:,1)=momentOut.varDecHorStates(:,:,1); 
tableStru.dataMat(:,:,2)=momentOut.varDecHorStates(:,:,4); 
tableStru.dataMat(:,:,3)=momentOut.varDecHorStates(:,:,8); 
tableStru.dataMat(:,:,4)=momentOut.varDecStates; 

tableStru.rowLabels=report.stateIRF;
tableStru.colLabels=shockNamesStru.short;
tableStru.corner={'States'};
latexStru.fileName=['Variance Decomposition States',strdate];
latexStru.maxCols =15;
latexStru.maxRows =20;
multiPageTable(tableStru,latexStru); 


%% 8.d Volatiliies
cd(cucd); 

% Volatility simulated observables
in.modenames = {''};
in.varnames  = Znames;
in.savename  = 'stdSimObs';
in.savepath  = outpath;
plot_stdhist(momentOut.STDSimObs, Y,in);

% Volatility simulated states
if ~isempty(report.stateMom)
    in.varnames = report.stateMom;
    in.savename = 'stdSimStates';
    plot_stdhist(momentOut.STDSimStates, [],in);
end

% Asymptotic Volatility of Observations
in.modenames = {''};
in.varnames  = Znames;
in.savename  = 'stdObs';
in.savepath  = outpath;
plot_stdhist(momentOut.STDObs, Y,in);

%% Asymptotic Correlations
% Observables
[ ac_data, ~, ~] = accov(Y,0:momentStru.NCorrel);
in.modenames = {'data','model'};
in.varnames = Znames;
in.savepath = outpath;
in.savename = 'crossCorrelationsObs';
plot_all_autocorr(momentOut.corrSimObs,ac_data,in);


% cd(outpath);
% save workspace;
% cd(cucd);




%%

[vdcom,vcvth,sdevf,ACzero,ACone,vcshocks,auone,vdech,vdech_st,vdcom_st]=vardecom_mode(GG,RR,SDX,ZZ,50);
var_outcell  = crtcell( vdcom, names.observables , shonames ,'Theoretical Variance Decomposition');
ac0th_outcell= crtcell( ACzero, names.observables, names.observables   ,'Theoretical AC(0)'); 
ac1th_outcell= crtcell( ACone , names.observables, names.observables   ,'Theoretical AC(1)'); 
cd(cucd); 
stdvsdat_outcell= crtcell([std(Y)' sdevf],names.observables  ,{'Data Std','Theoretical Std'},'Standard Deviation');  
ac1vsdat_outcell= crtcell([diag(accov(Y,1)) diag(ACone)],names.observables,{'Data AC(1)','Theoretical AC(1)'},'AC(1)');  

theoretical_cell=merge_cells(topcell,stdvsdat_outcell,2); 
theoretical_cell=merge_cells(theoretical_cell,ac1vsdat_outcell,2); 
theoretical_cell=merge_cells(theoretical_cell,var_outcell,2); 
theoretical_cell=merge_cells(theoretical_cell,ac0th_outcell,2); 
theoretical_cell=merge_cells(theoretical_cell,ac1th_outcell,2); 
clear var_outcell  ac0th_outcell ac1th_outcell stdvsdat_outcell ac1vsdat_outcell; 

par_cell = [parnames, num2cprec(parMode, 5)];
ss_cell = [ssnames, num2cprec(ssvec)];
mode_cell = merge_cells(par_cell, ss_cell, 2);
mode_cell = merge_cells(topcell, mode_cell, 2);
mode_cell = merge_cells(mode_cell, stdcell(:, 1:3), 2); 

%return 

%% 9. Save Moments of MODE 
cd(outpath)
if ~isunix
    xlswrite(mode_filename, mode_cell,'Details');
else
    savecell(mode_cell,[],outpath,['Details_', mode_filename]);
end

clear mode_cell; 
try
    if ~isunix
        xlswrite(mode_filename,spec_outcell,'spectrum');
    else
        savecell(spec_outcell,[],outpath,['spectrum','_', mode_filename]);
    end
    disp('Saved Spectrum'); 
    
catch
    disp('Spectrum Tab not saved')
end

try
    if ~isunix 
        xlswrite(mode_filename,theoretical_cell,'TheoreticalMoments');
    else
        savecell(theoretical_cell,[],outpath,['TheoreticalMoments', '_', mode_filename]);
    end
    
    disp('Save Theoretical_Cell'); 
    clear theoretical_cell; 
catch
    disp('Stationary Variance Tab not saved')
end

cd(cucd); 
flag_Tablesdone  =1; 
flag_Complete    =0; 

% cd(outpath); 
% save workspace; 
% cd(cucd); 
close all; 

estima_states; 


cd(outpath); 
save workspace; 
cd(cucd); 
disp('Hessian Done'); 

estima_out_analysis 
